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ABSTRACT 

Inactive  high  area-to-mass  ratio  (HAMR)  resident  space  objects  (RSOs)  near  the  geosynchronous  or¬ 
bit  (GEO)  regime  pose  a  hazard  to  active  GEO  RSOs.  The  combination  of  solar  radiation  pressure 
(SRP)  and  solar  and  lunar  gravitational  perturbations  cause  variations  in  the  orbit  parameters  of  the  in¬ 
active  RSOs.  The  high  A/m  nature  of  these  objects  results  in  greater  sensitivity  to  the  SRP  reflected  in 
mean  motion,  inclination  and  eccentricity.  The  subsequent  drift  with  respect  to  Earth-based  tracking 
sites,  combined  with  a  temporal  orientation  variability  with  respect  to  the  sun,  results  in  an  inability 
to  successfully  reacquire  many  of  these  RSOs  as  they  transition  through  periods  of  days  to  weeks  out 
of  view  of  observing  sites.  The  unknown  material  properties,  approximations  in  the  shadow  function 
used  for  eclipse  modeling,  errors  in  size  estimates  and  variations  in  the  orientation  with  respect  to 
the  sun  result  in  unpredictable  (random)  errors.  The  characteristics  of  the  error  distributions  are  of 
interest,  as  most  commonly  used  orbit  estimation  techniques  assume  a  Gaussian  distribution.  This 
work  examines  improvements  in  the  trajectory  estimation  when  applying  an  approach  that  attempts 
to  approximate  the  true  probability  density  function  (pdf)  of  the  state  error  distribution  instead  of 
assuming  it  to  be  Gaussian.  These  improvements  are  evaluated  in  the  context  of  operational  tracking 
capabilities. 

1.  BACKGROUND  AND  GOALS  OF  THE  STUDY 

Analysis  and  simulation  results  show  that  the  distributions  of  propagated  state  errors  resulting  from 
each  of  the  observation  and  modeling  uncertainties  to  be  non-Gaussian  [1].  Fig.  1.  shows  the  mean 
total  position  error  distribution  for  a  A/m  (AMR)  =  10  m2/kg  GEO  object  propagated  over  1-day, 
where  the  errors  include  unmodeled  rotation,  unmodeled  atmospheric  refraction  and  absorption  of  the 
solar  rays  during  a  1-hour  eclipse,  a  1%  error  in  the  estimated  AMR  and  1%  errors  in  the  specular 
and  diffuse  reflection  coefficients.  The  distributions  of  the  position  errors  from  the  individual  error 
contributions  are  shown  to  be  quite  different  from  the  total  position  error  distribution  only  due  to 
AMR  state  estimation  errors,  shown  in  Fig.  2.  In  summary,  the  results  of  the  study  show  that  the  SRP 
modeling  errors  result  in  the  non-Gaussian  distribution  of  position  and  velocity  errors. 

Results  of  recent  work  [2]  show  promise  when  an  adaptive  Gaussian  mixture  (AGM)  approach  is 
applied  to  the  determination  of  attitude  in  the  presence  of  attitude  observation  and  modeling  errors. 
The  purpose  of  this  work  is  to  take  an  initial  step  towards  understanding  the  utility  and  possible  ben¬ 
efits  of  implementing  an  AGM  approach  to  the  problem  of  trajectory  estimation  given  that  previous 
work  indicates  the  tendency  of  state  errors  becoming  non-Gaussian.  Since  the  AGM  methodology  is 
directly  aimed  at  approximating  the  true  state  distribution,  it  is  capable  of  yielding  a  more  realistic 
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Fig.  1.  Total  Position  Error  Distribution  for  a  HAMR  GEO  Object:  Combined  Rotation,  AMR, 
Specular  and  Diffuse  Reflection  Modeling  Errors 
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Fig.  2.  Total  Position  Error  Distribution  for  a  HAMR  GEO  Object:  AMR  Estimation  Errors  Only 

representation  of  the  state  distribution.  As  a  consequence,  the  more  realistic  state  error  uncertain¬ 
ties  captured  by  this  approach  should  translate  into  improved  accuracies  of  subsequent  HAMR  RSO 
reacquisition,  data/track  association,  and  conjunction  assessments. 

2.  SRP  MODELING  ERROR  INFLUENCE  ON  HAMR  OBJECT  ORBIT  ERRORS 

The  analysis  presented  in  this  paper  examines  the  state  estimation  errors  resulting  from  improper 
modeling  of  the  RSO  characteristics.  The  non-conservative  forces  affecting  RSOs  are  driven  by  these 
characteristics  and  given  that  HAMR  objects  have  unknown  characteristics  a  priori ,  the  research 
attempts  to  assess  operational  tracking  performance  in  the  presence  of  dynamics  mismodeling  effects. 
This  tracking  performance  assessment  is  done  by  comparing  a  method  that  is  constrained  to  Gaussian 
assumptions  and  one  that  is  not.  The  dynamic  models  focus  on  debris  objects  near  the  geosynchronous 
orbit  regime  having  AMR  values  in  the  range  of  0.1  —  20m2/kg.  The  gravitational  perturbations  are 
limited  to  the  dominant  zonal  harmonic  (J2)  and  luni-solar  perturbations  since  the  unknown  HAMR 
dynamic  effects  are  driven  solely  by  non-conservative  forces.  The  sensitivity  to  SRP  modeling  errors 
given  by  the  relatively  high  AMR  values  are  of  interest  in  assessing  the  orbital  prediction  sensitivity 
to  mismodeling  of  SRP  related  parameters. 


The  modeling  error  analyses  are  facilitated  by  comparing  trajectories  generated  via  a  “truth” 
model  which  uses  high  fidelity  force  models  versus  simpler  models  that  are  typically  used  in  the 
reduction  and  prediction  of  orbits  since  there  is  no  method  to  incorporate  a  physically  realistic  model 
with  a  lack  of  a  priori  information.  The  high  fidelity  models  implement  a  detailed  Earth  shadowing 
model  that  accounts  for  atmospheric  refraction  [3],  RSO  orientation  dependence,  and  thermal  emis¬ 
sions  acceleration  components  that  have  time- varying  material-dependent  thermal  properties.  Lastly, 
knowledge  of  the  exact  AMR  value  is  often  poorly  known,  and  so  errors  attributed  to  this  uncertainty 
are  modeled. 

The  detailed  SRP  acceleration  model  used  for  this  work  incorporates  both  reflective  and  absorptive 
effects  of  the  radiation  incident  on  the  RSO,  in  addition  to  illumination  with  respect  to  the  sun  which 
is  dependent  on  both  Earth  eclipsing  and  RSO  orientation  changes  of  the  modeled  surfaces.  The 
model  accommodates  any  number  of  specified  flat  surfaces  according  to  the  equations 

1  NS 

asRP  =  VI  Ami(Fp®Ri(ki  ■  hi)  +  Pt)  (1) 

c  i=i 

where  the  specular  and  diffuse  reflection  terms  for  each  surface  are  defined  as 

2 

Rt  =  (1  -  Si)ki  -  -dihi  +  2 (hi  •  (2) 

the  emission  term  is  defined  as 

Pi  =  ^dioTfhi  (3) 

The  parameter  c  is  the  speed  of  light,  a  is  the  Steffan-Boltzman  constant,  and  the  parameters  specific 
to  a  surface  i  are  the  specular  reflectivity  coefficient,  si5  the  diffuse  reflectivity  coefficient,  di,  the 
absorption  coefficient,  a*,  and  the  surface  temperature  T,.  The  incident  solar  flux  parameter,  <f>,  is 
adjusted  according  to  the  distance  of  the  space  object  from  the  sun,  but  has  a  nominal  (average)  value 
of  1367  Watts/m2  at  distance  of  1  Astronomical  Unit  (AU).  The  size  parameter  for  each  surface 
is  defined  as  the  area  to  mass  ratio  Am*  =  Arj rn  for  surface  area  A,  and  mass  m.  The  two  unit 
vectors  k,  and  n,  represent  the  incident  sun  and  surface  normal,  respectively,  for  surface  i.  Their 
dot  product  defines  the  attitude  dependence  of  the  SRP  model.  The  surface  normal  unit  vectors  are 
defined  in  a  body  reference  frame,  and  must  be  transformed  to  the  inertial  reference  frame.  Finally, 
the  parameter  Fp  (the  shadow  function),  in  the  interval  [0, 1],  is  considered  in  this  analysis  being  that 
it  is  the  function  used  to  model  passage  into  and  out  of  Earth  eclipse.  The  simplest  geometric  models 
are  that  of  a  cylinder  or  a  conic/fractional,  while  the  more  complex  pseudo-physical  and  physical 
models  attempt  to  account  for  the  effects  of  atmospheric  refraction  and  absorption  of  the  incident 
light  rays. 

The  temperature  variations  are  a  function  of  the  illumination  angle,  the  time  since  shadow  entry 
and  its  associated  illumination  angle,  and  the  time  since  shadow  exit  as  well  as  its  associated  illumi¬ 
nation  angle.  Here,  the  illumination  angle  is  defined  as  the  angle  between  the  incident  sun  vector  and 
a  vector  normal  (perpendicular)  to  a  given  surface.  From  the  radiation  pressure  modeling  work  done 
to  support  precision  orbit  determination  of  the  Topex/Poseidon  mission  [4]  the  truth  trajectory  models 
the  temperature  variations  resulting  from  the  illumination  transitions  according  to  assumed  materials 
properties. 

3.  ESTIMATION  STRATEGIES 
3.1  Problem  formulation 


To  best  understand  the  innovation  presented  in  this  work,  we  begin  with  the  general  formulation  of  the 
state  model  that  constitutes  the  foundation  of  the  orbit  determination  and  prediction  process.  Consider 


the  n-dimensional  continuous  dynamical  system  described  by  the  differential  equation 


x{t)  =  f{t,x{t))  +  G(t,x{t))w{t)  (4) 

where  x(t)  G  M"  is  the  state  vector,  /(•)  G  Mn  is  the  dynamical  system  model,  w(t)  G  Ms  is  the 
dynamical  system  process  noise,  and  G(-)  G  Mnxs  is  a  state-dependent  shape  matrix.  The  process 
noise  is  assumed  to  be  a  white-noise  process  with  mean  and  covariance 

E{iu(f)}  =  0  V  t  and  E  [w(t)wT(r)  j  =  Q(t)5(t  —  r) 

where  Q(t)  G  Msxs  is  the  process  noise  power  spectral  density  and  8(t  —  r)  is  the  Dirac  delta.  The 
continuous  dynamical  system  is  accompanied  by  m-dimensional  discrete-time  observations  described 

by 


zk  =  h(tk,  xk)  +  Lkvk  (5) 

where  zk  G  Mr"  is  the  measurement,  h  ( ■ )  G  M'"  is  the  observational  model,  vk  G  Rr  is  the  measure¬ 
ment  noise,  and  Lk  G  Mmx  r  is  a  shape  matrix.  The  measurement  noise  is  assumed  to  be  a  white-noise 
sequence  with  mean  and  covariance 

E  {^fc}  =  0  V  /c  and  E  {vkvk,}  =  Rk8k,k' 

Here,  5k,k'  represents  the  Kronecker  delta. 

Given  the  system  described  by  Eqs.  4  and  5,  the  goal  is  to  reconstruct  the  unknown  and  hidden 
dynamical  process,  x(t),  via  inference  utilizing  the  observed  measurements,  zk,  coupled  with  the 
assumed  dynamical/observational  relationships,  /(•)  and  h(-).  This  goal  must  be  accomplished  in 
the  presence  of  uncertain  dynamics  with  the  aid  of  simultaneous  noisy  measurements  from  multiple 
unknown  targets.  Therefore,  target  estimates  must  be  initiated  and  maintained  while  simultaneously 
associating  incoming  data  to  existing  estimated  targets  (or  classifying  the  data  as  belonging  to  new 
targets  in  order  to  initiate  new  estimates). 

In  the  context  of  the  problem  addressed  in  this  paper,  x(t)  is  the  6-element  vector  of  Cartesian  po¬ 
sition  and  velocity  to  be  estimated  in  the  orbit  determination  (OD)  process,  f(t.  x(t))  is  the  nonlinear 
model  of  the  forces  acting  on  the  object  (a  function  of  the  state),  and  w(t)  represents  the  modeling 
errors  present  in  the  nonlinear  forces.  Note,  the  process  noise  is  assumed  to  be  a  white  noise  process, 
though  in  reality  this  is  likely  not  the  case.  It  should  also  be  noted  that  the  SRP  coefficient  would  ide¬ 
ally  be  a  state  parameter  to  be  estimated.  Though  future  work  will  include  this,  the  SRP  model  error 
differences  incorporated  into  this  analysis  serve  to  address,  implicitly,  the  impact  of  SRP  estimation 
errors  and  their  associated  distributions. 

3.2  Track  initialization 

For  any  statistically  based  estimation  strategy,  a  method  for  obtaining  a  distribution  of  the  initial  state 
is  required.  One  approach  utilized  in  many  orbit  determination  applications  is  to  perform  an  ini¬ 
tial  orbit  determination  (IOD)  step  and  then  refine  the  initial  estimate  provided  by  IOD  via  a  Batch 
Processor  mechanism  (e.g.  a  square-root  information  filter  [5])  before  transitioning  to  sequential  es¬ 
timation  (e.g.  an  extended  Kalman  filter  or  an  unscented  Kalman  filter  [6,  7]).  The  classical/standard 
methods  for  performing  IOD  are:  Laplace’s  Method  [8],  Gauss’  Method  [8],  Double  r-Iteration  [8], 
or  Gooding’s  Method  [9].  All  of  the  aforementioned  IOD  methods  require  the  use  of  measurements 
comprised  of  six  independent  parameters  in  order  to  produce  a  six-parameter  orbit.  If  restrictions 
are  made  as  to  the  class  of  orbits  (i.e.  removing  some  of  the  parameters),  the  requisite  number  in¬ 
dependent  measurement  parameters  can  potentially  be  reduced.  For  example,  three  angle  pairs  (e.g. 
right  ascension  and  declination)  are  typically  required  for  determination  of  an  elliptical  (eccentric  ) 


orbit.  If,  however,  we  assume  the  orbit  is  circular,  only  two  pairs  (four  independent  measurements) 
are  required. 

The  approach  considered  here  is  based  on  the  notion  of  hypotheses  discussed  in  [10]  and  assumes 
an  angles-only  observational  scenario.  The  utilization  of  hypotheses  for  the  generation  of  state  esti¬ 
mates  allows  the  bypassing  of  the  IOD  state  and  facilitates  the  direct  entry  to  sequential  estimation. 
In  general,  if  a  measurement  of  the  right  ascension  and  declination  of  an  object  is  available  then  four 
more  quantities  are  required  to  create  a  state  vector  since  the  state  belongs  to  a  six-dimensional  space. 
However,  in  [10]  restricted  classes  of  orbits  are  considered,  allowing  for  the  use  of  fewer  than  four 
quantities  in  the  generation  of  hypotheses.  As  a  consequence  of  facilitating  a  more  broad  class  of 
orbits,  the  set  of  quantities  used  in  hypothesis  generation  for  this  work  is  taken  to  be:  the  geocentric 
range,  the  velocity  magnitude  along  the  radial  direction,  the  velocity  magnitude  along  the  tangential 
direction,  and  an  angle  defining  the  tangential  direction. 

Briefly,  the  conversion  from  measurements  and  hypotheses  to  a  state  is  presented.  The  right  ascen¬ 
sion,  a,  and  the  declination,  5,  are  made  available  via  an  optical  observation.  Given  a  hypothesized 
geocentric  range,  r,  the  inertial  position  of  the  object  with  respect  to  the  geocenter  is  constructed  by 
first  forming  the  unit  vectors 


Ustn  ~ 


stn 


and 


u 


stn  I 


obj / stn 


cos  5  cos  a 
cos  S  sin  a 
sin  5 


which  represent  the  unit  vector  of  the  observing  site  with  respect  to  the  geocenter  and  the  unit  vector 
of  the  object  with  respect  to  the  observing  site.  The  angles  7  and  (3  are  then  computed  via 

(  1 1  ^  1 1 

7  =  vr  -  cos ~l{ulobj/stn  ■  u\tn )  and  (3  =  ir  -  7  -  sin^1  (  — ^ -  sin  7 

which  facilitates  the  direct  computation  of  the  range  from  the  observing  site  to  the  object  as 


P 


2r||r*J|  cos  (3 


The  hypothesized  position  of  the  object  with  respect  to  the  geocenter  may  then  be  found  to  be 


obj  P'U'obs/stn  ^ stn  (^) 

Now,  using  the  hypothesized  object  position  and  the  hypothesized  angle  defining  the  tangential  direc¬ 
tion  (i.e.  6),  the  unit  vector  pointing  to  the  object  with  respect  to  the  geocenter  and  a  rotation  matrix 
are  computed  as 

ur  =  ‘fJ  and  T  —  J3x3  —  sin0[wrx]  +  (1  —  cos#)[ttrx]2 

WKbjW 

where  I3x3  represents  the  identity  matrix  of  square  dimension  three  and  [-x]  represents  the  skew- 
symmetric  matrix  resulting  from  the  vector.  T  is  a  rotation  by  angle  6  about  axis  ur,  which  represents 
the  orientation  of  the  tangential  (in-track)  direction  and  is  bounded  by  the  orbit  inclination.  The 
hypothesized  in-track  unit  vector  is  then  given  by 


Ui  =  T 


e3  x  u, 

1 63  X  ur 


where  e3  e  M3  has  the  third  element  to  be  unity  and  the  remaining  elements  to  be  zero.  In  the  case 
where  bme3  =  ur,  another  unit  vector  can  be  used  in  place  of  e3  without  loss  of  generality.  Finally, 


the  hypothesized  velocity  of  the  object  with  respect  to  the  geocenter  may  be  computed  making  use  of 
the  hypothesized  radial  velocity,  vr,  and  the  hypothesized  tangential  velocity,  vg,  as  follows 

Vlobj  =  VrUr  +  VgUi  (7) 

Therefore,  given  the  measurement  data  as  a  and  <5  and  the  hypothesis  data  as  r,  6,  vr,  and  vg,  a 
hypothesized  position  and  velocity,  i.e.  Eqs.  6  and  7,  of  the  object  with  respect  to  the  geocenter  may 
be  determined.  This  process  is  repeated  for  each  combination  of  hypotheses. 

3.3  Multiple  hypothesis  filter 

The  multiple  hypothesis  filter  (MHF)  is  essentially  a  bank  of  UKF  filters  which  are  weighted  based 
on  the  time-history  statistical  agreement  to  incoming  data.  The  initial  bank  of  filters  is  initiated 
by  transforming  each  hypothesis  into  a  state  mean  and  covariance  via  the  unscented  transform  [7] 
which  allows  the  measurement  uncertainty  and  uncertainty  associated  with  the  hypothesis  selection 
to  be  transformed  into  state  uncertainty.  The  filters  are  initialized  with  equal  weights.  As  more  data 
are  processed,  the  MHF  recursively  refines  the  available  hypotheses,  updating  the  filter  weights  and 
pruning  hypotheses  which  have  not  agreed  to  acquired  data.  When  one  hypothesis  dominates,  the 
MHF  becomes  equivalent  to  a  single  UKF,  making  it  a  winner-take-all  method.  However,  it  is  also 
possible  that  multiple  hypotheses  are  consistent  with  the  collected  data  due  to  the  existence  of  multiple 
local  minima  or  a  geometry-induced  lack  of  observability,  or  that  more  than  one  hypothesis  converges 
to  the  same  local  (or  global)  minimum.  More  specific  information  on  the  MHF  methodology  and  its 
implementation  can  be  found  in  [10]. 

3.4  Adaptive  Gaussian  mixtures 

At  the  end  of  a  measurement  track,  the  MHF  algorithm  is  typically  left  with  one  or  two  Gaussian 
hypotheses  that  propagate  to  the  next  measurement  using  an  unscented  Kalman  filter  to  maintain  the 
mean  and  covariance.  The  AGM  methodology  attempts  to  characterize  the  time  evolution  of  a  state 
probability  density  function  (pdf),  rather  than  a  simple  mean  and  covariance,  by  assuming  that  the 
distribution  can  be  represented  by  a  finite  sum  of  Gaussian  components.  These  components  are  cho¬ 
sen  based  on  the  hypotheses  after  they  are  conditioned  by  measurements.  For  stochastic  continuous 
dynamic  systems,  the  exact  evolution  of  the  state  pdf  is  given  by  the  Fokker-Planck- Kolmogorov 
(FPK)  equation  [11,  12], 

which  is  a  formidable  problem  to  solve  because  of  the  following  issues: 

-  Positivity  of  the  pdf:  p(t,  x)dx  >  0. 

-  Normalization  constraint  of  the  pdf:  jRn  p(t,  x)dx  =  1. 

-  No  fixed  Solution  Domain-,  how  to  impose  boundary  conditions  in  a  finite  region  and  restrict 
numerical  computations. 

Note  that  in  Eq.  8  and  in  the  following  developments,  the  derivative  of  a  scalar  with  respect  to  a 
vector  has  been  taken  to  be  a  column  vector  and  not  a  row  vector.  Analytical  solutions  to  the  FPK 
equation  are  restricted  to  a  limited  class  of  dynamic  systems,  and  numerical  methods  suffer  in  higher 
dimension  systems.  Approximating  the  forecast  state  pdf  using  a  finite  Gaussian  mixture  model  offers 
a  way  to  impose  the  aforementioned  restrictions  on  the  solution  of  the  FPK  equation,  and  leads  to  a 
convex  quadratic  programming  problem  for  which  straightforward  solutions  exist. 


Consider  the  following  equation  depicting  the  Gaussian  mixture  model  approximation  for  the 
forecast  conditional  density  function,  p(t,  x(t)  j  Z\:) 

N 

P(t,  x(t)  |  Zk)  =  J]  wlt]kpi  with  Pi  =  A f(x(t)  |  \x\ k,  PI ]fc)  (9) 

1=1 


Here,  and  PE  represent  the  conditional  mean  and  covariance  of  the  ith  component  of  the  Gaus¬ 
sian  mixture  pdf  with  respect  to  the  k  measurements,  and  w^k  denotes  the  amplitude  of  ith  Gaussian  in 
the  mixture.  Furthermore,  J\f  (x  \  fi,  P)  represents  the  calculation  of  a  Gaussian  probability  density 
function  for  the  random  variable  x  with  mean  p,  and  covariance  P,  i.e. 


N  (x  \  ix,  P ) 


1 


exp 


p)TP~1(x 


The  positivity  and  normalization  constraints  on  the  mixture  pdf,  p(t,  x(t)  \  Zk),  lead  to  constraints  on 
the  amplitudes  as 


N 

X/aV  =  1  and  WA  k  >  0  v  t  (10) 

i= 1 

In  [13],  it  is  shown  that  since  all  the  components  of  the  mixture  pdf  of  Eq.  9  are  Gaussian,  only 
estimates  of  their  mean  and  covariance  need  to  be  maintained  between  tk  and  tk+l  in  order  to  obtain 
the  optimal  state  estimates.  These  can  be  propagated  by  using  a  continuous-time  derivation  of  the 
unscented  Kalman  filter  [6,  14],  which  uses  a  set  of  deterministically  chosen  sigma-points  that  capture 
the  mean  and  covariance  of  the  initial  distribution. 

Notice  that  the  weights  Wi  corresponding  to  each  Gaussian  component  are  also  unknown.  Hence, 
the  idea  is  to  use  FPK  equation  error  as  a  feedback  to  update  the  amplitude  of  different  Gaussian 
components  in  the  mixture  pdf.  In  other  words,  we  seek  to  minimize  the  FPK  equation  error  under 
the  assumption  of  Eqs.  9  and  10.  Substituting  Eq.  9  in  Eq.  8  leads  to 

e{t,x{t))  =  X ^  -  Zk'>  CFPK(p{t,x{t)  |  Zk))  (11) 

where  Cfpk(-)  is  the  FPK  operator.  The  application  of  this  operator  is  given  by, 

N 

£ FPK(p(t,x(t )  |  Zk))  =  Wl\k^FPK(Pi )  (12) 

i= 1 


whith  the  application  of  the  FPK  operator  to  the  individual  Gaussian  components  given  by 

=  /(«,*(*)) (13) 

+  Itrac  e{QM9x(t)^(i)T} 

The  first  term  in  Eq.  1 1  is  the  time  derivative  of  the  probability  density  function  approximated  by  the 
Gaussian  mixture,  and  is  given  by, 


dp{t,x{t)  |  Zk) 


N 

E 

2—1 


+  w 


t\k 


dpi 

dp, 


l 

t\k 


T 

tft\k  +  ^Vtrace 


dt 


(14) 


The  derivatives  of  the  first  two  moments  are  given  by  UKF  equations  [14],  and  the  derivative  of  the 
weights,  w>i,  is  obtained  by  the  following  finite  difference  approach: 


<\k  =  ^  «| k  ~  wl\k )  where  t'  =  t  + At  (15) 

One  point  of  interest  is  the  computation  of  tangent  linear  dynamics  (the  dynamics  Jacobian)  re¬ 
quired  for  use  in  Eq.  13,  i.e.  the  computation  of 


F(t,x(t)) 


df(t,x(t)) 

dx(t) 


Since  the  UKF  is  being  utilized,  the  tangent  linear  dynamics  matrix  is  not  directly  available  and  a 
method  for  approximating  this  matrix  must  be  employed.  Consider  the  cross-covariance  between  the 
state  and  the  dynamics  of  the  state,  represented  by  Pxx.  In  the  case  of  linear  dynamics  with  respect  to 
the  state,  the  cross-covariance  may  be  computed  as 


Pxx  =  PxFT(t,x(t )) 


which  implies  that  if  the  covariance  on  the  state  is  known  and  the  cross-covariance  between  the  state 
and  the  dynamics  is  known,  then  the  tangent  linear  dynamics  may  be  approximated  as 

F(t,x(t))  =  JVP-1 

Utilizing  the  UKF,  the  cross-covariance  is  found  to  be 

2  n 

Pxx  =  y^w\c)(Xj  -  mx)(Xi  -  mx)T 
i= 0 

(c) 

where  n  is  the  dimension  of  the  state,  w\  1  are  covariance  weights  [14],  X,  are  the  sigma-points, 
mx  is  the  mean  value  of  x ,  X ,  are  the  transformed  sigma-points  given  by  evaluating  the  dynamical 
system  at  each  of  the  sigma-points,  and  mx  is  the  mean  value  of  the  transformed  sigma-points.  Thus, 
implementing  a  UKF  structure,  the  tangent  linear  dynamics  matrix  can  be  approximated  an  applied 
to  Eq.  13. 

Now,  at  a  given  instant  in  time  after  propagating  the  mean  and  covariance  of  each  Gaussian  com¬ 
ponent,  we  seek  to  update  the  weights  of  the  Gaussian  mixture  model  by  minimizing  Eq.  1 1  over  the 
volume  of  interest.  This  requires  the  computation  of  several  integrals  that  lead  to  a  convex  quadratic 
programming  problem  with  a  unique  solution,  of  which  the  full  derivation  can  be  found  in  [15].  Us¬ 
ing  this  formulation,  one  can  approximate  the  time-evolution  of  the  conditional  state-pdf  at  any  time 
instant. 

A  major  challenge  in  solving  the  convex  minimization  problem,  however,  is  the  need  to  evaluate 
integrals  involving  Gaussian  pdfs  over  volume  V .  This  volume  V  can  be  defined  making  use  of  the 
approximation  that  the  mass  of  a  Gaussian  pdf  is  concentrated  in  a  finite  volume  around  its  mean. 
This  is  one  of  the  major  advantages  of  using  Gaussian  mixture  models,  because  the  space  over  which 
probability  mass  lies  is  easy  to  define.  The  full  derivation  of  these  integrals  can  be  found  in  [15],  and 
they  can  be  computed  exactly  for  polynomial  nonlinearity.  In  general,  numerical  integration  methods 
are  necessary  such  as  Gaussian  quadrature,  Monte  Carlo  integration  or  unscented  transformation  [16]. 
While  in  lower  dimensions  the  unscented  transformation  and  Gaussian  quadrature  methods  are  mostly 
equivalent,  the  unscented  transformation  is  computationally  more  appealing  in  higher  dimensions 
because  the  number  of  points  taken  to  evaluate  the  integral  grows  only  linearly  with  the  number  of 
dimensions. 


3.5  Probabilistic  data  association 


In  situations  involving  simultaneous  measurements  of  multiple  targets,  association  of  the  incoming 
data  to  existing  tracks  is  an  integral  step  in  computing  accurate  solutions.  There  are  many  data  asso¬ 
ciation  schemes  discussed  in  the  literature;  however,  the  method  used  in  this  work  is  the  probabilistic 
data  association  method  (PDA). 

The  method  employs  a  statistical  generalization  of  the  typical  Euclidean  distance,  termed  the 
Mahalanobis  distance,  which  is  computed  by 

d2  =  (zk  -  z~)Tp-\zk  -  z~) 


where  zk  is  the  incoming  measurement  data,  is  the  predicted  value  of  the  measurement  data,  and 
Pz  is  the  covariance  of  the  difference  between  the  incoming  and  predicted  measurement  data  (i.e.  the 
innovations  covariance).  It  can  be  shown  that  for  Gaussian  distributions  the  Mahalanobis  distance 
satisfies  a  y2  distribution.  Specification  of  a  probability  gate  can  then  be  used  to  validate/associate 
(or  can  also  be  thought  of  as  a  scheme  to  reject)  measurements  with  a  y2  distribution  table.  However, 
if  several  measurements  are  validated/associated,  then  there  exists  a  hard  decision.  That  is,  it  must  be 
decided  which  measurement(s)  to  associate  and  which  measurement(s)  to  reject. 

The  PDA  method  [17]  offers  a  data  association  solution  which  avoids  the  need  to  make  a  hard 
decision  and  instead  fuses  all  relevant  information  based  on  statistical  agreement.  Given  m  available 
measurements  at  time  tk,  the  measurements  are  determined  to  be  valid  if  their  computed  Mahalanobis 
distance  falls  below  some  threshold  dictated  by  a  probability  gate  [17].  The  number  of  validated 
measurements,  v,  is  such  that  v  <  m.  For  each  of  the  associated  measurements,  the  association  event 
probabilities  are  calculated  as 


<%•  = 


U(> 


e;,w(4i  %-,&) 


V 

such  that  ^  a.i  =  1 

i— 1 


The  association  event  probabilities  describe  how  likely  it  is  that  a  validated  measurement  actually 
originated  from  the  target  under  consideration.  Once  the  association  event  probabilities  have  been 
determined,  each  of  the  validated  measurements  may  be  processed  to  produce  an  updated  distribution, 
which  is  computed  as  a  mixture  of  the  individual  updates  from  each  validated  measurement  [10,  17]. 


4.  SIMULATION  AND  ANALYSIS  RESULTS 
4.1  Simulation  description 

The  overall  simulation  environment  consists  of  modules  which  initialize  HAMR  objects  according 
to  input  distribution  characteristics  for  the  orbit,  object  shape,  and  material  properties,  integrate  the 
HAMR  objects  over  a  specified  period  of  time,  simulate  measurements  of  right-ascension  and  dec¬ 
lination,  and  then  apply  estimation  strategies  to  attempt  to  recover  the  object’s  orbit.  To  integrate 
the  HAMR  objects  over  a  period  of  time,  a  high-fidelity  model  accounting  for  the  accelerations  due 
to  gravity  (point-mass,  zonal  harmonics,  and/or  spherical  harmonics),  third  body  perturbations  ac¬ 
counting  for  accurate  positioning  of  the  bodies  (especially  solar  and  lunar  perturbation),  and  solar 
radiation  pressure  with  attitude-dependent  and  thermal  radiation  effects  is  applied.  In  generating  the 
simulated  measurements,  methods  for  the  accounting  for  observer  lighting  conditions,  target  object 
lighting  conditions,  light-time  transit  delay  [18],  stellar  aberration  effects  [19],  and  sensor  field-of- 
view  limitations  are  employed.  Furthermore,  a  specified  distribution  for  a  measurement  bias  and  a 
noise  are  sampled  to  add  corruption  to  the  measurements.  In  addition,  a  Monte  Carlo  capability  is  im¬ 
plemented  to  generate  different  measurement  sequences  by  sampling  the  measurement  bias  and  noise 
distributions  on  each  run.  This  capability  provides  a  mechanism  for  generating  filtering  solutions  for 
many  different  measurement  sequences  as  a  way  of  testing  robustness  to  the  measurements. 


The  measurement  availability  for  the  scenario  considered  is  summarized  by:  an  initial  tracking 
segment  of  15  minutes  with  measurements  every  30  seconds,  a  data  gap  lasting  6  hours,  a  final  track¬ 
ing  segment  of  15  minutes  with  measurements  every  30  seconds,  and  a  final  data  gap  lasting  17.5 
hours.  The  measurements  utilized  are  topocentric  right  ascension  and  declination  from  a  station  with 
a  known  location  with  each  angle  measurement  having  a  zero-mean  random  constant  bias  with  a 
standard  deviation  of  0.5  arc-seconds  and  a  zero-mean  random  noise  with  a  standard  deviation  of  1 
arc-second. 

4.2  Orbit  and  parameter  distributions 

In  the  analysis,  a  truth  reference  orbit  is  generated  via  numerical  integration  with  specified  orbit 
parameter  (Keplerian  element)  distributions  for  the  initial  conditions  of  the  integrator  and  by  speci¬ 
fication  of  SRP-related  material  properties.  The  Keplerian  element  distributions  considered  for  the 
analysis  are  all  uniform  distributions  with  specified  minimum/maximum  parameter  values,  which  are 
taken  to  be 

38,  000  [km]  <  a  <  43,  000  [km] ,  0  <  e  <  0.4 ,  0  [deg]  <  %  <  10  [deg] 

0  [deg]  <  0  <  360  [deg] ,  0  [deg]  <  c o  <  360  [deg] ,  and  0  [deg]  <  M  <  360  [deg] 

where  a  is  the  semi-major  axis,  e  is  the  eccentricity,  i  is  the  inclination,  Q  is  the  right-ascension  of  the 
ascending  node,  to  is  the  argument  of  perigee,  and  M  is  the  mean  anomaly.  The  material  properties  are 
generated  by  randomly  selecting  a  material  with  a  60%  probability  of  the  material  being  end-of-life 
(EOL)  MLI  Kapton,  a  30%  probability  of  the  material  being  solar  cells,  and  a  10%  probability  of  the 
material  being  EOL  white  paint.  The  selection  of  the  material  then  determines  the  diffuse  reflectance, 
d,  and  the  specular  reflectance,  s.  The  values  used  in  this  analysis  are  [20] 


EOL  MLI  kapton  : 

d  = 

0.03 

and 

s  =  0.35 

solar  cells  : 

d  = 

0.04 

and 

s  =  0.04 

EOL  white  paint  : 

d  = 

0.53 

and 

s  =  0.03 

In  all  cases,  the  absorption,  a,  is  chosen  so  that  the  three  material  property  values  sum  to  unity,  i.e. 
a  +  d  +  s  —  1.  Finally,  the  area-to-mass  ratio  is  selected  randomly  using  a  uniform  distribution  with 
specified  minimum/maximum  value,  which  are  taken  to  be 

0.1  [m2/kg]  <  —  <  20  [m2/kg] 

TTL 

4.3  MHF  and  AGM  results  comparison 

A  Monte  Carlo  analysis  sampling  the  measurement  bias  and  noise  distributions  on  each  run  was 
performed  for  the  MHF  algorithm.  The  estimation  errors  and  the  associated  3cr  estimation  error  co- 
variances  are  shown  in  radial,  in-track,  and  cross-track  (RIC)  coordinates  in  Fig.  3.  and  Fig.  4.  In 
both  figures,  the  consistent  behavior  of  the  MHF  is  clearly  demonstrated.  The  algorithm  is  consis¬ 
tently  able  to  initiate  the  state  estimate,  refine  that  initial  estimate,  reacquire  the  object  after  a  period 
of  no  data,  and  then  forward  predict  the  state  estimates.  However,  it  is  also  seen  that  the  collection 
of  Monte  Carlo  runs  is  not  very  representative  of  the  filter’s  perception  of  the  estimation  errors  with 
clearly  present  biases  and  non-normal  behavior.  Furthermore,  it  can  be  seen  in  Fig.  4.  that  at  some 
point  during  the  simulation  run,  the  estimation  error  from  every  Monte  Carlo  run  strays  beyond  the 
3<t  covariance  bound  of  the  filter,  indicating  that  the  filter  is  not  properly  capturing  the  behavior  of 
the  estimation  error. 

Since  the  PDA  method  for  data  association  relies  on  an  accurate  accounting  of  the  state  estimation 
errors  in  order  to  properly  gate  and  accept/reject  measurement  data,  it  is  possible  that  the  MHF  will 


Fig.  3.  Estimation  error  of  Monte  Carlo  runs  for  the  MHF  algorithm  for  the  full  24  hours  of  the 
simulation.  Estimation  error  is  in  blue  and  the  3 a  estimation  error  covariance  is  shown  in  red. 
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Fig.  4.  Estimation  error  of  Monte  Carlo  runs  for  the  MHF  algorithm  for  the  time  period  starting 
after  the  end  of  the  second  measurement  arc.  Estimation  error  is  in  blue  and  the  3cr  estimation  error 
covariance  is  shown  in  red. 


fall  into  a  situation  in  which  measurements  become  consistently  rejected.  This  can  be  clearly  seen 
from  Fig.  4.  at  14  hours  past  epoch  in  the  cross-track  component  where  all  of  the  samples  are  outside 
of  the  3a  covariance  bound.  An  attempt  to  reacquire  the  object  at  this  time  would  result  in  a  failure 
of  the  algorithm  to  properly  associate  data  due  to  the  misrepresentation  of  the  state  uncertainty  by  the 
MHF.  To  help  solve  this  problem,  an  AGM  was  implemented  in  an  attempt  to  better  characterize  the 
state  estimation  errors.  The  MHF  and  AGM  methods  were  applied  to  a  single-run  case  in  order  to 
demonstrate  the  relative  performance  of  the  two  algorithms.  The  RIC  state  estimation  errors  and  the 
associated  filter  covariances  for  the  two  methods  are  shown  in  Fig.  5.  and  Fig.  6.  It  is  seen,  especially 


in  Fig.  6.  that  the  state  estimation  errors  for  the  two  methods  are  comparable;  however,  the  major 
difference  lies  in  the  difference  between  each  filter’s  perception  of  the  potential  distribution  of  the 
errors.  That  is,  the  AGM  3a  covariance  is  sufficient  enough  to  encompass  the  estimation  error  at  the 
points  where  the  MHF  3a  covariance  does  not  encompass  the  MHF  estimation  error.  As  previously 
mentioned,  this  is  a  particularly  important  aspect  in  the  association  of  measurement  data  to  existing 
estimates.  With  a  more  proper  representation  of  the  state  estimation  error  distribution  data  association 
can  function  more  appropriately,  which  would  lead  to  a  better  reacquisition  of  objects. 
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Fig.  5.  Estimation  error  for  the  MHF  and  AGM  algorithms  for  the  full  24  hours  of  the  simulation. 
Estimation  error  for  the  AGM  is  in  solid  blue  and  for  the  MHF  is  in  dashed  green.  The  3cr  estimation 
error  covariance  for  the  AGM  is  in  solid  red  and  for  the  MHF  is  in  dashed  black. 


Finally,  to  compare  the  measurement  prediction  at  the  end  of  the  time  span,  contours  of  constant 
probability  were  computed  for  both  the  MHF  and  AGM  method.  These  are  plotted  for  two  different 
probability  values  in  Fig.  7.  Along  with  the  contours,  the  true  measurement  as  well  as  the  prediction 
of  the  measurement  value  for  the  MHF  and  the  AGM  are  plotted  for  comparison.  The  curves  of  con¬ 
stant  probability  are  computed  by  taking  the  peak  value  of  the  pdf  for  each  method  and  determining 
the  locus  of  points  corresponding  to  a  scaled  value  of  this  peak  probability,  i.e.  50%  of  the  peak  in 
Fig.  7(a).  and  1%  of  the  peak  in  Fig.  7(b).  The  percentage  of  peak  probability  corresponds  to  con¬ 
fidence  intervals  with  the  smaller  percentage  corresponding  to  a  larger  confidence  region.  That  is,  in 
defining  an  area  upon  which  to  search  for  the  object,  the  area  defined  by  a  50%  peak  is  going  to  be 
smaller  than  the  area  defined  by  a  1%  peak.  From  Fig.  7(a).,  it  can  be  seen  that  for  the  50%  contour, 
the  MHF  contour  does  not  encompass  the  true  measurement  while  the  AGM  contour  does  encompass 
the  true  measurement.  In  the  case  of  a  1%  contour,  both  methods  have  contours  encompassing  the 
true  measurement.  In  both  figures,  it  is  clearly  seen  that  the  AGM  prediction  of  the  measurement  is 
closer  than  that  of  the  MHF.  Ultimately,  the  fact  that  a  contour  corresponding  to  a  larger  percentage 
of  peak  probability  encompasses  the  result  means  that  a  smaller  relative  area  can  be  searched  in  an 
attempt  to  locate  a  previously  tracked  HAMR  RSO. 

5.  CONCLUSIONS  AND  FUTURE  WORK 

Monte  Carlo  analysis  of  an  MHF  method  applied  to  the  initialization,  tracking,  and  prediction  of  a 
HAMR  object  has  been  illustrated  and  shown  to  exhibit  non-Gaussian  behaviors  and  the  potential 
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Fig.  6.  Estimation  error  for  the  MHF  and  AGM  algorithms  for  the  time  period  starting  after  the  end 
of  the  second  measurement  arc.  Estimation  error  for  the  AGM  is  in  solid  blue  and  for  the  MHF  is  in 
dashed  green.  The  3 a  estimation  error  covariance  for  the  AGM  is  in  solid  red  and  for  the  MHF  is  in 
dashed  black. 


(a)  50%  of  peak  probability 


(b)  1%  of  peak  probability 


Fig.  7.  Final  measurement  pdf  contours.  Note:  axes  are  scaled  differently  to  show  detail. 


for  loss  of  tracking.  The  implementation  of  an  AGM  method  to  the  HAMR  object  has  also  been 
illustrated,  showing  promise  in  its  ability  to  more  realistically  represent  the  state  estimation  error 
knowledge  (covariance).  Results  presented  here  have  shown  the  following: 

-  The  state  error  distributions  for  HAMR  RSOs  are  non-Gaussian. 

-  A  Monte  Carlo  analysis  shows  how  the  state  estimation  errors  exceed  the  reported  uncertainties 
from  the  estimation  strategy.  This  clearly  shows  that  the  covariances  reported  are  optimistic 
and  if  they  were  to  be  used  for  data/track  association,  could  lead  to  data/track  rejections  (i.e. 
false  negatives).  If  there  were  closely  space  objects,  it  could  lead  to  false  positives  (i.e.  the 
data/track  association  associates  two  separate  objects  believing  them  to  be  the  same  one). 


-  The  AGM  confirms  a  non-Gaussian  state  estimation  error  pdf  given  two  significantly  weighted 
components. 

-  The  AGM  uncertainties  always  encompass  the  state  errors  (more  accurate/realistic  than  MHF 
alone). 

-  The  AGM  predicted  observations  are  closer  to  the  actual  observation  than  those  predicted  by 
the  MHF  alone  (more  accurate  mean  than  MHF  alone). 

These  results  indicate  the  potential  for  allowing  reacquisition  methods  to  function  more  properly 
due  to  the  more  realistic  characterization  of  the  uncertainty.  As  data  outages  become  longer,  the  need 
for  accurate  characterizations  of  the  uncertainties  becomes  more  pressing  for  the  success  of  followup 
data  association  and  tracking. 

Future  work  will  focus  on  utilizing  the  MHF  and  AGM  techniques  for  identifying  the  necessary 
requirements  on  the  angles-only  data  quality  in  order  to  guarantee  data/track  association  and  reacqui¬ 
sition  of  HAMR  RSOs.  Specific  hypotheses  to  be  addressed  are  as  follows: 

-  An  MHF  augmented  by  an  AGM  is  able  to  better  characterize  RSOs  given  poor  or  unknown  a 
priori  information  than  traditional  estimation  methods. 

-  Using  an  AGM  can  significantly  improve  HAMR  RSO  reacquisition  over  traditional  estimation 
strategies,  especially  in  the  presence  of  sparse  angles-only  data. 

-  For  HAMR  RSO  data/track  association  and  characterization,  the  quantity  of  angles-only  data 
can  be  significantly  reduced  if  fused  with  temporal  photometric  signatures 

-  The  AGM  is  able  to  significantly  aid  in  HAMR  RSO  ID/discrimination  as  compared  to  tradi¬ 
tional  estimation  strategies. 

-  The  AGM  does  is  better  at  approximating  the  true  pdf  as  compared  to  traditional  estimation 
strategies. 

-  The  more  realistic  state  estimation  error  distributions  produced  by  AGM  improves  the  realism 
of  conjunction  assessments. 

In  order  to  guarantee  data/track  association  and  reacquisition  of  HAMR  RSOs,  the  angles-only 
data  quantity,  duration  and  quality  required  must  be  analyzed,  with  sensor  noise  and  biases  limited 
to  known  performance  values.  The  ultimate  goal  of  new  processes  that  provide  improvements  to 
data/track  association  and  conjunction  assessment  will  first  be  demonstrated  via  simulation  analysis, 
and  eventually,  using  actual  observational  data  either  serendipitously  collected,  or  obtained  through 
coordinated  observational  experiments. 
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